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ABSTRACT 

We use a set of high-resolution hydrodynamical simulations of clusters of galaxies to 
study the build-up of the intracluster light (ICL), an interesting and likely significant 
component of their total stellar mass. Our sample of groups and clusters includes 
AGN feedback and is of high enough resolution to accurately resolve galaxy popula- 
tions down to the smallest galaxies that are expected to significantly contribute to the 
stellar mass budget. We describe and test four different methods to identify the ICL in 
cluster simulations, thereby allowing us to assess the reliability of the measurements. 
For all of the methods, we consistently find a very significant ICL stellar fraction 
(~ 45%) which exceeds the values typically inferred from observations. However, we 
show that this result is robust with respect to numerical resolution and integration 
accuracy, remarkably insensitive to changes in the star formation model, and almost in- 
dependent of halo mass. It is also almost invariant when black hole growth is included, 
even though AGN feedback successfully prevents excessive overcooling in clusters and 
leads to a drastically improved agreement of the simulated cluster galaxy population 
with observations. In particular, the luminosities of central cluster galaxies and the 
ages of their stellar populations are much more realistic when including AGN. In the 
light of these findings, it appears challenging to construct a simulation model that 
simultaneously matches the cluster galaxy population and at the same time produces 
a low ICL component. We find that intracluster stars are preferentially stripped in 
a cluster's densest region from massive galaxies that fall into the forming cluster at 
z > 1. Surprisingly, some of the intracluster stars also form in the intracluster medium 
inside cold gas clouds that are stripped out of inf ailing galaxies. 

Key words: galaxies: clusters: general - galaxies: formation - cosmology: theory - 
methods: numerical - black hole physics 



1 INTRODUCTION 

In optical observations, galaxy clusters appear as concen- 
trations of galaxies on the sky. However, deep exposures 
reveal that not all the emission comes from stars that re- 
side in the cluster's member galaxies. Instead, there is also 
a smoothly distributed stellar component which is typically 
peaked around the cluster's central galaxy, but extends to 
large radii. Due to its low surface brightness, observations of 
this intracluster light (ICL) component are difficult, result- 
ing in a significant uncertainty in the current observational 
constraints of the amount of ICL present in clusters. 

iLin and Mohrl \2004 ) give an overview of the fractions 
of intracluster stars reported in the literature for different 
clusters and groups (see their Table 2). The reported values 
span a huge range from ~ 2% to ~ 50%. Inste ad of analyz- 
ing ve ry deep exposures of individual objects, IZibetti et al.l 
l|2005h stacked hundreds of cluster images from the Sloan 
Digital Sky Survey and masked all satellite galaxies, allow- 



ing them to obtain a very deep exposure of the brightest 
cluster galaxy (BCG) and the ICL of an average galaxy 
cluster. For a mean cluster mass of 7 — 8 x 1O 13 M0 they 
find that about ~ 22% of the stars reside in the BCG, while 
roughly ~ 1 1% are intraclus t er stars with in a fixed aperture 
of 500 kpc. iGonzalez et all (|2005l . 120071 1 describe another 
study for which a significant number of clusters and groups 
were observed. By analyzing the observed systems individu- 
ally, they found that both the total stellar fractions as well 
as the relative fractions of stars residing in the BCG+ICL 
component strongly increase with decreasing halo mass. It 
was also concluded that intracluster stars are important for 
the baryon budget of clusters and groups. 

The latter is particularly interesting in the light of 
the recent findings that the baryon fractions measured in 
clusters seem to be significantly lower than those inferred 
from the most curre nt cosmological constraints (see e.g. 
iMcCarthv et al.ll2007l h However, no mechanism that can ef- 
ficiently segregate baryons from dark matter on the scale 
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of massive galaxy clusters is known. On the other hand, if 
there is indeed a significant component of intracluster stars, 
correctly accounting for them may relax the reported ten- 
sions. 

The distribution and abundance of intracluster stars 
were also investigated based on cosmological hydrody- 
namical simulations of galaxy cluster formation that in- 
cluded radiative cooling and heating, as well as models 
for st a r formation and supe rn ova feedback llMurante et al.l 
2004 IWillman et al.1 |2004 ISommer-Larsen et alj 1200*1 



Murante et all 120071 : iDolag et all I2009I ). The origin of 
intracluster stars w as studied in particular detail in 
iMurante et al.1 l|2007t ) finding an ICL fraction that rises with 
halo mass from ~ 10% at the group scale to ~ 30% for mas- 
sive clusters. In their simulations, about half of the intraclus- 
ter stars came from galaxies associated with the merger tree 
of the BCG, and most intracluster stars were liberated from 
their former host galaxies during merger events after redshift 
z = 1. Using sim ilar simulations bu t a different method to 
identify the ICL, IDolag et al.l (|2009l ) obtained average ICL 
fractions of around ~ 33%, with a significant scatter on the 
group scale. No trend of the mean ICL fraction with halo 
mass was detected when their new ICL definition was used. 

A major problem of most hydrodynamical cluster sim- 
ulations thus far has been that they suffered from exces- 
sive "overcooling" within the densest cluster regions, where 
the gas cooling times are short. As a result, an unrealis- 
tically large fraction of cold gas and consequently a large 
amount of stars in clusters formed from the strong cool- 
ing flows. This in turn leads to central galaxies that are 
too bright and too blue, as well as extremely high total 
stellar fractions. It is widely believed that some source of 
non-gravitational energy input, like heating by active galac- 
tic nuclei (AGN), is necessary to offset the cooling in cluster 
cores and t o make simulated clu sters compatible with obser- 
vations. In lSiiacki et al.l (|2007h and iPuchwein et al.l (2008) 
it was shown that including a model for AGN feedback in 
such simulations, indeed, strongly improves the agreement 
between simulated a nd observed X-ray prop erties of clusters 
and groups (see also lMcCarthv et aljBoOSh . 

In this work, we explore whether AGN feedback can 
also resolve discrepancies between simulated and observed 
cluster galaxy populations, and in particular, between the 
amount of intracluster light. Our high-resolution simulation 
set is well suited to investigate the origin of the ICL because 
for the first time a large simulated sample of galaxy clusters 
and group is available which can resolve the cluster galaxy 
populations down to very low galaxy masses and at the same 
time includes a successful model for the growth and feedback 
activity of AGN. 

This paper is structured as follows. In Section 2 we de- 
scribe our simulation set, the numerical models used, and 
the analysis methods we employ for measuring the ICL com- 
ponent in the simulations. We also summarize the tests we 
have carried out to investigate the robustness of our results 
with respect to numerical parameters and variations in the 
star formation model. In Section 3, we present our main 
results on the origin and amount of ICL found in our simu- 
lations, as well as on the properties of the simulated cluster 
galaxy populations. Finally, we summarize our findings and 
give our conclusions in Section 4. 



2 METHODS 

We have carried out high-resolution cosmological hydrody- 
namical simulations of a large sample of galaxy clusters and 
groups, including a treatment of star formation and feed- 
back processes. We will here mainly focus on analyzing the 
distribution of stars in the formed objects. To this end we 
develop and intercompare several distinct methods to mea- 
sure the following components in our simulated clusters: the 
central galaxy (BCG)lj the satellite galaxies (i.e. all other 
cluster member galaxies), and the ICL. We then investigate 
the properties of these components and their relation to each 
other. 



2.1 The simulations 

Our simula ted galaxy cluster and group sample is described 
in detail in lPuchwein et ail l^OOcf ). In this work we use the 
16 most massive objects of the sample. They approximately 
uniformly cover the mass range from M^jjf = 2xl0 13 M© to 
1.5xlO 15 M 

In brief, we have s elected dark matter h alos from the 
Millennium simulation (|Springel et alj|2005bh at z = and 
resimulated them at higher mass and force resolution, in- 
cluding a gaseous component and accounting for hydrody- 
namics, radiative cooling, heating by a UV background, star 
formation and supernovae feedback. For each halo, at least 
two kinds of resimulations were performed. One containing 
the physics just described, and an additional one that also 
included a model for feedback from active galactic nuclei 
(AGN). We note that the original halo selection was only 
based on mass and otherwise random. New initial condi- 
tions for the resimulations were created by populating the 
Lagrangian region of each halo in the original initial condi- 
tions with more particles and adding additional small-scale 
power, as appropriate. At the same time, the resolution has 
been progressively reduced in regions that are sufficiently 
distant from the forming halo. Gas has been introduced into 
the high-resolution region by splitting each parent particle 
into a gas and a dark matter particle. 

We have adopted the same flat ACDM cosmology as 
in the parent Millennium simulation, namely fi m = 0.25, 
Qa = 0.75, h — 0.73, n s — 1 and erg = 0.9. A baryon density 
of Ob = 0.04136 has been chosen for consistency with the 
cosmic bary on fraction inferred from current cosmological 
constraints dKomatsu et alJl200Sl ). 

The simulations were run with the GADGET-3 code 



1 The central galaxy, i.e. the galaxy closest to the minimum of 
the cluster potential well, is typically also the brightest cluster 
galaxy, except for very rare exceptions. For simplicity, we always 
refer to the central galaxy of a simulated cluster or group when 
we use the abbreviation BCG. 

2 Throughout this paper, spherical overdensity masses and radii 
will be denoted by the symbols M and r, where we use the su- 
perscripts "crit" or "mean" to indicate whether the overdensity 
is with respect to the critical density of the Universe or the mean 
matter density of the Universe at the cluster redshift, respectively. 
The subscripts 200 and 500 indicate how many times larger the 
mean density inside the corresponding spherical region is com- 
pared to the chosen reference density. 
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zoom 


softening 


particle mass [h 




factor 


[h^kpc] 


DM 


gas 


stars 


1 


7.5 


8.3 x 10 s 


1.8 x 10 s 


8.9 x 10 7 


2 


3.75 


1.1 x 10 8 


2.1 x 10 7 


1.0 x 10 7 


3 


2.5 


3.1 x 10 7 


6.2 x 10 6 


3.1 x 10 6 


4 


1.875 


1.3 x 10 7 


2.8 x 10 6 


1.4 x 10 6 



Table 1. Force and mass resolutions used in the cluster and 
group simulations. The physical (Plummer-equivalent) softening 
lengths that were used for redshifts z < 5 are listed in Column 
2. For z > 5, comoving softenings of A5h~ kpc, 22.5/i _1 kpc, 
15 h~ 1 hpc, and 11.25 /i —1 kpc were used for zoom factors 1. 2, 3, 
and 4, respectively. The masses of the dark matter (DM), gas and 
star particles in the simulations are tabulated in Columns 3-5. 



(based on ISpringej 120051 ). which employs an entropy- 
conserving formulation of smoothed particle hydrodynam- 
ics (SPH). Radiative cooling and heating was calculated 
for an optically thin plasma of hydrogen and helium, and 
for a time-varying but spatially uniform UV background. 
Star formation and supernovae feedback were modelled 
with a subreso lution multi-phase mode l for th e interstellar 
medium as in ISpringel and Hernquistl |2003al ). Black hole 
gr owth and associated f eedb ack processes were f ollowed as 
in lSpringel et al.l (|2005al ) and lSiiacki et al.l (|2007h . 

Table [T] summarizes the different mass and force res- 
olutions that were used in our resimulations. For all halos 
with Mfoo ^ 1-5 x 10 14 h~ 1 MQ simulations with a zoom fac- 
tor of 3 (see Table [TJ have been performed, while all more 
massive cluster were resimulated only with a zoom factor 
of 2. For some objects several resimulations have been per- 
formed, covering a wide range in mass and force resolution. 
This allows us to assess the convergence of the simulations. 



2.2 Identifying intracluster stars and cluster 
galaxies 

Starting from the particle distribution in the simulated clus- 
ters, we aim to identify the cluster's central galaxy, all satel- 
lite member galaxies and the ICL component. As a first 
step, we run a friends-of-friends (FoF) group finder with 
a linking length of 0.2 in units of the mean interparti- 
cle distance, as well as the SUBFIND substructure finder 
(|Springel et al.l [200ll V The FoF algorithm is applied only 
to the high-resolution dark matter particles, with the gas 
and stars linked to their nearest dark matter particle. This 
avoids biases in the group catalogue due to the very differ- 
ent spatial distribution of dark matter and baryonic parti- 
cles. We then use a version of SUBFIND that was specifically 
adapted to work r obustly with hydrodynamical simulations 
l|Dolag et al.l 120081 ) for decomposing each FoF-group into a 
main halo and self-bound substructures. We consider the 
stellar content of each self-bound substructure of the cluster 
to be a satellite galaxy. All star particles that are part of 
the main halo are assumed to belong either to the BCG or 
to the ICL component. It is, however, not straightforward 
to unambiguously decide to which of these two components 
such a star particle should be assigned. In the following, we 
therefore employ several independent methods for making 



a distinction between the BCG and the ICL in simulated 
clusters. These are: 

• Method 1: A cut-off radius around the central galaxy. 
The simplest possible approach is to just cut off the BCG 
at a specific radius r cut . However, the size of BCGs corre- 
lates with cluster mass. We therefore choose to scale r cut 
with cluster mass. We obtain an appropriate scaling re- 
lation by combining an empiri cal relation between clus- 
ter mass and BCG luminosity (jPopesso et al.l [2007) with 
a relation between B CG luminosity and half-light radius r c 
l|Bernardi et al.ll2007h . We decide to cut at r cut = 3 x r c , 
yielding 

/ A/crit \ 0.29 

When using this method, all star particles of the main halo 
that are within the cut-off radius are considered to be part 
of the BCG, while all main halo star particles outside this 
radius are considered to be intracluster stars. 

• Method 2: Analysis o f the surface brightne ss profiles. 
A de Vaucouleurs profile l|de Vaucouleurs! Il948) is usually 
a good fit to the surface brightness profile of an elliptical 
galaxy. For cD galaxies, however, one often finds a light ex- 
cess at large radii with respect to such a fit to the inner re- 
gion l|Schombertlll98rJ| . This light excess can be attributed 
to intracluster stars. 

The complete surface brightness profile of BCG and ICL 
can typically be well fit by the sum o f two de Vaucouleurs 
profiles (see e.g. [Gonzalez et al.l [20051 ). We perform such a 
two-component fit to the surface brightness profile of all 
stars residing in a simulated cluster's main halo, as found 
by SUBFIND, and consider the component with the smaller 
effective radius as the BCG and the component with the 
larger effective radius as the ICL. This procedure is illus- 
trated for a cluster with M^o = 2.2 x 10 14 h' 1 M Q in 
Fig. [1] We note that we obtain projected mass profiles for 
the BCG and the intracluster stars by assuming that in each 
radial bin both components have the same mass-to-light ra- 
tio. In order to compute surface brightness profiles we use 
the stellar population synt hesis model library GALAXEV 
IIBruzual and Chariot! [2003) to assign a r-band luminosity 
to each star particle according to its mass, age and metal- 
licity. 

• Method 3: Analysts of the stellar velocity distribution. 
In relaxed simulated clusters, the distribution of the ve- 
locities v of the main halo stars can be well fit by a sum 
of tw o Maxwell-Boltzmann distributions (see Dola g et al.l 
120091) . Here v is given by v — u s tar - u c ius|, where i> s tar 
is the velocity of the star particle and i> c i us is the veloc- 
ity of the centre of mass of the cluster's main halo. In this 
method, we perform such a fit for each cluster and define 
the BCG mass as the mass of the stars in the Maxwell- 
Boltzmann component with the lower velocity dispersion, 
while we associate the higher dispersion component with 
the ICL. This is illustrated in Fig. [2]for a cluster with mass 
M^oo = 2.2 x 10 14 h~ 1 M Q . 

• Method 4 : Analysis of the stellar velocity distribution 
and binding energies. While Method 3 allows an estimate of 
the total stellar BCG mass, it does not assign individual star 
particles to either the BCG or ICL component, as it yields, 
for each velocity bin, only the fraction of stars in each com- 
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ponent. However, by also including the binding energies of 
the star particles in the analysis, the method can be ex- 
tended to allow such an ind i vidua l assignment, as discussed 
in full detail in lDolag etafl l|2009l ). 

In short, we use this approach in the following way. We 
start with the two component Maxwell-Boltzmann fit ob- 
tained from Method 3. Then we use an iterative algorithm 
that assigns star particles according to their binding energy 
either to the BCG or to the ICL, with the iteration repeated 
until the BCG and ICL velocity distributions roughly match 
the fitted Maxwell-Boltzmann components. We start the it- 
eration with some value r grav for the size of the BCG's grav- 
itating mass distribution. Then, the gravitational potential 
of all simulation particles with a distance r < r grav from 
the cluster centre is computed. All main halo star particles 
that are bound in this potential are assigned to the BCG in 
this iteration step, while unbound main halo stars are as- 
signed to the ICL. Next, the velocity dispersion of the ICL 
component is calculated and compared to the velocity dis- 
persion of the two Maxwell-Boltzmann fits. If the velocity 
dispersion of the ICL is too large, 7- g r av is reduced to un- 
bind more of the slowly moving particles. Otherwise, r grav 
is increased. We stop the iteration once the ICL velocity dis- 
persion agrees with the velocity dispersion of the higher dis- 
persion Maxwell-Boltzmann component to within 1%. Fig. [2] 
displays the resulting velocity distributions of the BCG and 
ICL components obtained in this way and compares them 
to the corresponding Maxwell-Boltzmann fits. 

In Fig. [T] we compare the dynamical definition of BCG 
and ICL obtained with Method 4 to the analysis of the sur- 
face brightness profiles used in Method 2. For the shown 
cluster the agreement between these two completely different 
methods is remarkably good. However, in some other clus- 
ters the differences are admittedly larger (see also Fig. (6|. 

We note that we can obtain a BCG and ICL mass by all 
of these four methods. However, only Method 1 and Method 
4 allow us to individually assign each main halo star particle 
to one of these two components. 

2.3 Robustness with respect to numerical 
resolution and integration accuracy 

Before we apply the methods introduced above to our full 
cluster and group sample, we investigate whether our satel- 
lite galaxy, BCG and ICL stellar masses are robust with re- 
spect to numerical resolution by analysing three halos that 
were simulated at several different mass and force resolu- 
tions. We also test our integration accuracy with one fur- 
ther cluster simulation where we used two times smaller time 
steps and a higher force accuracy setting in our tree code(f| 
Fig. shows the masses of all stars, of all main halo 
stars (i.e. BCG+ICL), and of the BCG in these test simu- 
lations for a group with Mago = 4.4 x 10 13 h~ 1 M Q and a 
cluster with mass M^oo = !- 5 x 1014 h~ 1 M Q . We computed 
these masses within r§oo, based on resimulations with zoom 
factors 1, 2, and 3, without inclusion of AGN physics. The 
BCG masses obtained with all the four different methods 
introduced in Sect. [231 are shown. The total stellar masses 

3 In this run, the opening angle for nodes in the tree part of the 
force computation is a factor ~ 1.6 smaller. 




1 2 5 10 20 50 100 200 500 
radius r [kpc/h] 

Figure 1. r-band surface brightness profiles of main halo stars of 
a 2.2 X 10 14 7i -1 M e cluster. The red solid curves show the surface 
brightness profile of all main halo stars and its two-component de 
Vaucouleurs fit (smooth curve). The straight green dashed and 
blue dotted lines indicate the two individual de Vaucouleurs com- 
ponents corresponding to the BCG and ICL. The additional green 
dashed and blue dotted curves indicate the surface brightness pro- 
files of the BCG and ICL as found by Method 4- For this cluster 
the results of Method 2 and Method 4 agree very well. The r-axis 
has been chosen to be linear in r 1 ' 4 . The simulation has been 
performed without AGN feedback. 



are well converged even in our lowest resolution simulations 
without AGN. However, the mass of stars in satellite galax- 
ies, which is given by the difference between the results 
shown for all stars and for the BCG+ICL component, is 
not fully converged at the resolution of zoom factor 1. On 
the other hand, the results for zoom factors 2 and 3 for the 
BCG+ICL and the satellite galaxy mass agree very well, in- 
dicating that for zoom factor 2 and larger good convergence 
is achieved. 

The BCG masses that are inferred with the different 
analysis methods can sometimes slightly differ from each 
other. This will be discussed in more detail in Section [3+2] 
Here, we want to assess the numerical convergence of the re- 
sults obtained by each individual method. The BCG masses 
obtained by using a cutoff radius [Method 1) or the surface 
brightness profile analysis [Method 2) are already roughly 
converged even at the zoom factor 1 resolution. The results 
for Method 3 and Method 4 on the other hand, which both 
rely on fitting the velocity distribution of main halo stars, 
seem to be somewhat more easily affected by low numerical 
resolution, because the masses obtained for zoom factor 1 
are larger than those obtained from the higher resolution 
runs. However, at least for the group shown in the left panel 
of Fig. [3] this seems partly due to a slightly different timing 
of a rather massive subhalo that is just passing by close to 
the cluster center in the low resolution simulation, so it is 
not really clear how generic this effect is. In any case, we 
find that the BCG masses have robustly converged at zoom 
factors 2 and larger. 

However, we also find that the resolution requirements 
to obtain converged results for the total stellar mass in runs 
with AGN feedback are somewhat higher. This is because 
at higher resolution star formation shifts to smaller halos, 
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Figure 2. Velocity distribution of main halo stars of a 2.2 X 
10 14 h _1 Mq cluster. The red solid histogram shows the velocity 
distribution of all main halo stars. The smooth red solid curve 
is the two component Maxwell-Boltzmann fit used in Method 3 
and 4- The smooth green dashed and blue dotted curves indicate 
the individual Maxwell-Boltzmann distributions corresponding to 
the BCG and ICL. The green dashed and blue dotted histograms 
show the velocity distributions of the BCG and ICL as identified 
by Method 4- The simulation has been performed without AGN 
feedback. 
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Figure 3. BCG, BCG+ICL and total stellar mass found within 
r| U Q in a 4.4 X 10 13 /i _1 A/q group (left panel) and a 1.5 X 
10 14 ^"'Mq cluster (right panel) simulated without AGN at dif- 
ferent numerical resolutions. Results are shown for zoom factor 1, 
2, and 3 resolutions (see Table [TJ as well as for the four different 
methods for determining a BCG mass (see Sect. \2.2\ . In the right 
panel, two zoom factor 2 simulations are shown, one with the 
standard parameters (a) and one with smaller time steps and a 
more accurate force computation (b). Results are well converged 
for zoom factors 2 and larger. 



which contain smaller black holes that are less effective in 
suppressing star formation by their feedback. As a conse- 
quence, the total stellar mass in our simulated clusters is 
not nearly as well converged at the resolution of zoom fac- 
tor 1 as in our runs without AGN. It is, however, reasonably 
converged at zoom factor 2 (stellar mass about ~ 20% low) 
and well converged at zoom factor 3 resolution (about ~ 6% 
low). 

The results of the run with smaller time steps and 
higher force accuracy agree extremely well with those ob- 
tained form the corresponding run with our standard nu- 
merical parameters (see right panel of Fig. [3}. We hence 
conclude that the orbit integration is sufficiently accurate 
in our simulations. This is important and reassuring, as an 
insufficient integration accuracy may lead to an artificially 
enhanced stripping rate of stars out of cluster galaxies. 

In Fig. [4] we explore the convergence of the properties 
of the satellite galaxy population in more detail. Here the 
combined stellar mass of all satellite galaxies whose individ- 
ual stellar masses exceed M is shown as a function of M for 
simulations at different resolution. Comparing the curves ob- 
tained from simulations with zoom factors 1, 2, 3, and 4, we 
conclude that the satellite galaxy population is reasonably 
converged for stellar masses larger than ~ 5 X 10 10 h~ x MQ, 
~ 6 x 10 9 h~ 1 M Q and ~ 2 x 10 9 h _1 M Q , at zoom factors 1, 
2 and 3, respectively. In other words, the formation, merger 
and tidal disruption rates of galaxies above these masses 
should be robust with respect to numerical resolution. Us- 
ing the double Sche chter fit to the obse rved galaxy stellar 
mass function from Baldrv et al.l (|2008l ), we estimate that 
~ 69%, ~ 17% and ~ 9% of the total stellar mass are con- 
tained in galaxies with individual stellar masses below these 
resolution limits. Thus, at zoom factor 1 the majority of 
stars resides in poorly resolved galaxies. This explains why 
we obtain lower values for the mass in satellite galaxies in 
these runs (see Fig.[3J). On the other hand, in the runs with 
zoom factors 2 and 3, the vast majority of stars resides in 
galaxies whose abundance is numerically converged. We thus 
do not expect a significant spurious contribution to the ICL 
component form underresolved galaxies in these simulations. 

Overall, these test simulations indicate that the distri- 
bution of stars in our cluster resimulations with zoom factors 
2 and 3 are robust with respect to mass and force resolu- 
tion, except for maybe a slight underestimation of the total 
star formation efficiency in zoom factor 2 runs that include 
AGN. Importantly, we have also verified that our results are 
invariant when the time step size is reduced or the accuracy 
of the gravitational force computation is improved. 

2.4 Dependence on the star formation model 

We also investigated how sensitive our results are to changes 
in the parameters of the star formation model. More pre- 
cisely, we performed two test simulations, one with a higher 
gas density threshold for star formation, and another one in 
which we assumed a larger timescale for the conversion of 
gas into stars. 

The idea behind using a higher gas density threshold 
was that this might lead to smaller and denser galaxies, 
which should then be less prone to tidal stripping. This in 
turn could reduce the amount of ICL in the simulated clus- 
ters. In one of our tests, we therefore increased the thresh- 
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Figure 4. Combined stellar mass of all satellite galaxies with 
individual stellar masses exceeding M as a function of galaxy 
stellar mass M. Results are shown for zoom factor 1, 2 and 3 
1.5 X 10 14 h~ Mpi cluster and zoom 



simulations of a iW|no 
factor 3 and 4 simulations of a M"qq 



via 



3 X 10 13 h~ 1 M( 3 group. No 
AGN have been included in these runs. The vertical lines, arrows 
and, labels (zfl), (zf2), and (zfS) indicate the galaxy mass ranges 
in which the stellar mass in satellite galaxies is approximately 
converged in zoom factor 1, 2, and 3 simulations, respectively. 



old density for star formation by a factor of ~ 4 compared 
to the original simulation . Within the subresolution model 
of ISpringel and Hernquistl (|2003al ) for the regulation of star 
formation this was achieved by increasing the 'supernova 
temperature' Tsn from 10 8 K to 4 x 10 8 K, as well as the 
supernova evaporation factor Ao from 1000 to 4000, while 
keeping the gas consumption timescale the same. However, 
the analysis of this simulation yielded very similar results 
to the original run. The fraction of stars in the BCG+ICL 
component changed by only about 1%. 

Varying instead the timescale for star formation and 
thereby shifting the peak of the star formation activity to a 
different stage in the assembly history of the cluster might 
also affect the ICL component. In particular, since merg- 
ers play an imp ortant role in libera ting stars from their 
host galaxies fsee lMurante et al"1l2007l ), the efficiency of this 
mechanism for creating ICL stars may depend on the star 
formation history of the cluster. The idea behind our sec- 
ond test of the star formation model was therefore to in- 
crease the timescale to for the conversion of gas into stars 
from 2.1 Gyrs to 8.4 Gyrs. Keeping everything else equal, 
this required again to set Tsn = 4 x 10 8 K and Ao = 4000 
to leave the gas density threshold at its original value, while 
t$ itself was increased by a factor of 4. Using this larger gas 
consumption timescale indeed shifts the peak of the star for- 
mation history to much later times and broadens it at the 
same time. It also results in more star formation at low red- 
shift. However, also for this simulation we find remarkably 
little difference in the galaxy population and the ICL com- 
pared to the original run. Despite the substantial change 
in the star formation history, the fraction of stars in the 
BCG+ICL component increased by only 4%. 

Overall, these tests show that the fraction of stars that 
get unbound from their host galaxies and become intraclus- 



ter stars is remarkably insensitive to the detailed parameters 
of the star formation model. 



3 RESULTS 

We now turn to the presentation of our primary results, 
obtained by applying the methods introduced in Sect. 12.2 
to our full sample of simulated clustersQ In particular, wc 
investigate what fraction of a cluster's stellar mass resides 
in the BCG, the ICL, and in the satellite galaxies. We also 
study whether the answer to this question depends on cluster 
mass. Finally, we examine where the stars that end up in 
these three components at redshift z — are formed, and in 
what objects they fall into the cluster. 



3.1 Intracluster light, satellite galaxies, and BCGs 
in simulations and observations 

3.1.1 The baryonic mass fractions 

Figure [5] illustrates the baryon and stellar mass fractions of 
our simulated clusters, as well as the mass fraction of stars 
that reside in cluster galaxies as a function of cluster mass 
Mfoc i • Results from simulations both with and without AGN 
feedback are shown, a nd for comparis o n wit h observational 
constr aints , data from [Gonzalez et al.l (|2007l ) . rGiodini et al.l 
(|2009h and bin et ail (|2003h are included. 

In the simulations without AGN feedback, the baryon 
mass fractions within rfao differ by less than ~ 10% from the 
assumed cosmic baryon fraction. When AGN are included, 
this holds only for the most massive clusters, while there is 
a significant baryon depletion in poor clusters and g roups 
due to the AGN heating (see also iPuchwein et al.ll2008T ). 

Without AGN, we obtain very large stellar mass frac- 
tions. In massive clusters, roughly 40% of the baryonic mass 
is found in stars, while on the group scale this number rises 
to about 60%. This seems to be inconsistent with obser- 
vations, especially for massive clusters. When AGN are in- 
cluded, the stellar mass fractions within r§oo are lowered by 
about one third. One can compa re these total stellar m ass 
fractions to the constraints from iGonzalez etHl l|2007l) . as 
both cluster galaxies and intracluster stars were accounted 
for in that study. The simulated and observed stellar mass 
fractions are in good agreement on the group scale. However, 
we do not find a comparably s trong trend with halo mass 
in our simulations as inferred bv lGonzalez et all l|2007l '). and 
accordingly we obtain significantly larger stellar mass frac- 
tions for massive clusters. 

Also shown in Fig. [S] is the mass fraction of stars in 
cluster galaxies. More precisely, these values were calculated 
by summing up the stellar masses of the BCG, which was 
found by Method 1, i.e. using the radial cut-off scaled with 
the cluster mass, and of all satellite galaxies within r^oo- 
The resulting values can be compared to the mass fr a ctions 
of st ars in clus t er gal axies obtained bv lGiodini et al.l (|2009h 
and iLin et all (|2003h . Without AGN, we again find a sig- 
nificant discrepancy for massive clusters. Compared to the 
observations the combined stellar mass of the cluster galax- 
ies is too large. Also on the group scale the mass fractions in 



4 Unless indicated otherwise, we use the redshift 



outputs. 
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Figure 5. Mass fraction of baryons (circles), all stars (including 
ICL, squares), and stars residing in galaxies (triangles) within 
r 500 f° r clusters and groups simulated with and without AGN 
feedback. The latter was found by summing up the stellar masses 
of the BCG, which was found by Method 1, and all satellite galax- 



The I Gonzalez et al.l l|2007f ) data should be compared to the mass 
fraction of all stars, while the other observational data should 
be compared to the mass fraction of stars found in galaxies. The 
cosmic baryon fraction assumed in the simulations is shown for 
reference. The orange symbols show results of a simulation where 
also kinetic supernova feedback was included. 



the simulations are larger than those from lLin et al. I (l2003h. 
Howev er, they are roughly consistent with the Giodini et al.l 
(2009) results there. When AGN feedback is included the 
combined stellar mass of the cluster galaxies drops by about 
~ 40%. This substantially improves the agreement with ob- 
servations for massive clusters. On the g r oup s cale, there is 
then good agreemen t with the lLin et al. (|2003t ) data, while 
iGiodini et all (j2009T l inferred slightly larger mass fractions 
there. 

Overall, the stellar mass fractions obtained without 
AGN feedback are clearly too large. When the effects of 
AGN heating are included, the stellar mass in cluster galax- 
ies becomes however almost consistent with observations. 
In particular, t he total stellar mass, i .e. including the ICL, 
agrees with the iGonzalez et all (120070 results on the group 
scale. However, for massive clusters, we find a much larger 
ICL contribution in the simulations than they infer from ob- 
servations. The simulations also do not fully reproduce the 
observed strong trend of the stellar mass fractions with halo 
mass. 

One should mention that we did not include kinetic 
supernova feedback in the present set of simulations, in- 
stead all feedback from star formation was injected ther- 
mally. Using supernova feedback that includes a kinetic 



component can significantly reduce the amount of stars 
formed in hydrodynamical c o smological simulati ons (e.g 
ISpringel and HernquistJ l2003bl ; iBorgani et all 12006m . While 
this might help to reproduce the overall stellar mass frac- 
tions in m assive clusters, it is, considering that such simu- 
lations (see lBorgani et al]|2006t ) also find only a weak trend 
of stellar mass fractions with halo mass, unlikely to signif- 
icantly improve the overall agreement with observations 
Nevertheless, we repeated the simulation of one of our clus- 
ters, this time including a strong kinetic supernova feedback 
component in order to assess how it changes the fraction of 
stars found in the BCG, the ICL, and in the satellite galax- 
ies. This run also included AGN feedback. Its results are 
indicated by orange symbols in Fig. [5] The corresponding 
simulation without winds, is the one indicated by the red 
symbols at almost the same (slightly lower) M|oq. 



3.1.2 The stellar fractions in the BCG, ICL, and satellite 
galaxies 

In the following, we discuss what fraction of stars in our sim- 
ulated clusters and groups reside in the BCG, the ICL, and 
in the satellite galaxies, and how this compares to observa- 
tions. Fig. [6] shows the mass fractions of the stars found in 
the BCG and in the whole main halo, i.e. in the BCG+ICL 
component. The values were calculated within rfoo. Results 
obtained from simulations both with and without AGN are 
shown. The BCG fractions are reported for all the four 
different methods we used to determine a BCG mass (see 
Sect. 12 . 2 p . For reference, observatio nal constraints on the 
BCG+ICL luminosity fraction from IGonzalez et all (|20071 ) 
are also shown. They should be compared to the simulated 
BCG+ICL luminosity fractions indicated in the plot. The 
latter were computed within a projected radius (rather than 
a 3D radius) of r§oo , using the same method to obtain stellar 
luminosities as described in Sect. 12.21 

Without AGN, we obtain BCG+ICL fractions of 
roughly ~ 50% for massive clusters, and of about ~ 70% 
for groups. Note, however, that these values are sensitive to 
the radius inside which they are measured, because the dis- 
tribution of the BCG+ICL component is more concentrated 
than the distribution of the satellite galaxies. For example, 
within r^oo the values drop to about ~ 45% and ~ 65%, 
respectively. Accordingly, the BCG+ICL fractions are also 
lower when calculated within a projected radius rather than 
a 3D radius equal to r§oo- However, the BCG+ICL lumi- 
nosity fractions indicated in the figure are also affected by 
small differences in the mass-to-light ratio. Namely, we find 
mass-to-light ratios for the BCG+ICL components that are 
a few percent smaller (up to 10% smaller for massive clus- 
ter, while being almost identical for groups) than those of 
the whole stellar population in the simulated halo. Thus, the 



5 Recently, more sophisticated kinetic supernova feedback 
schemes that tune the wind velocity and mass loading fac- 
tor to properties of the host galaxy have bee n suggested 
llOppenheimer and Dave|[2006l ; lOkamoto et al.| [2009V Such mod- 
els might be more successful in reproducing a strong trend of the 
stellar mass fraction with halo mass. 

6 We assumed a supernova wind velocity of ~ 480 km/s and a 
wind mass flux rate that equa ls twice the star formation rate (for 
details of the wind model see[gpringel and Hcrnquist 20 03ah 
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Figure 6. Fraction of the stellar mass within r^g found in the 
BCG, and in the BCG+ICL component. Mass fractions are shown 
as a function of cluster mass M£qq and for the four different meth- 
ods of determining a BCG mass (see Sect. [2721) . The upper panel 
shows the values computed from runs without AGN, while the 
lower panel illustrates the results obtained when including the 
AGN feedback. Also shown is the BCG+ICL i-band luminosity 
fraction within a projected radius on the sky (rather than 3D ra- 
dius) of fgQQ. This data can be direct ly compared to the obs erved 
BCG+ICL luminosity fractions from Gonzale z et al. | 11200711 . 



i-band luminosity fractions are only slightly smaller than the 
mass fractions calculated within a 3D radius equal to r^oo- 
When accounting for these effects, we nevertheless still ob- 
tain larger B CG+ICL fractions th an derived from the ob- 
servations of iGonzalez et al] (|2007l ). especially for massive 
clusters. 

When AGN feedback is included, the stellar mass of the 
BCG+ICL component decreases by about ~ 25%. However, 
the total stellar masses decrease somewhat more strongly 
(see Fig. [5j). Thus, we find even slightly larger BCG+ICL 
fractions in the simulations with AGN. A potential explana- 
tion for this effect is that galaxies falling into the cluster in 
simulations with AGN may be more prone to tidal disrup- 
tion due to their lower stellar mass and due to an expansion 
of their halos as a consequence of the removal of gas by AGN 
heating. 

Especially in the runs without AGN, there are several 
clusters for which all the four different analysis methods 
for identifying the BCG yield very similar BCG fractions 
and masses (see also Fig. [1]). There are, on the other hand, 
also clusters where the results differ. In most of them, how- 
ever, there is still good agreement between the two methods 



that were inspired by observational approaches, i.e. between 
Method 1, where a radial cut-off is used, and Method 2, where 
surface brightness profiles are analyzed. Also, the methods 
that rely on the dynamics, i.e. Method 3 and Method 4, usu- 
ally agree well with each other, even though this is perhaps 
not too surprising as they both are based on the same fit of 
the stellar velocity distribution. In the rare cases were they 
do disagree, it is typically in unrelaxed objects, e.g. when 
an infalling group perturbs the stellar velocity distribution 
and prevents an accurate two-component Maxwell fit. 

Overall, there is hence some uncertainty in the BCG 
mass due to the choice of the analysis method used for sep- 
arating the BCG from the ICL. However, for all our differ- 
ent methods only a small part of the stars in the main halo, 
i.e. in the BCG+ICL component, may be assigned to the 
BCG. This also means that independent of how exactly we 
choose to make a distinction between BCG and ICL, we find 
a very significant ICL component in our simulated clusters 
and gro ups. We also note tha t the BCG and ICL compo- 
nents in IGonzalez etafl (|200rJ) are defined very similarly to 
our Method 2, hence our analysis method should yield results 
that can be directly compared to observations. 

The fraction of stars in the ICL can be read off form 
the difference between the BCG+ICL and BCG fractions in 
Fig. [5] As can be seen, the ICL fractions in our simulations 
are almost independent of the cluster mass. Within rjjoo, 
they are roughly ~ 45% without AGN and ~ 50% with 
AGN physics. If we go to the larger radius of rfoo instead, 
they are slightly lower, that is about ~ 40% without AGN 
and ~ 45% with AGN. 

At a first glance these findings s eem to be at odds with 
the results of iMurante et al.l (|2007l ) , where an ICL frac- 
tion rising with cluster mass was found, from ~ 15% at 
10 14 /i _1 M Q to ~ 30% at 1O 15 /i _1 M . However, the lower 
values they obtained are most likely a consequence of the 
much lower numerical resolution in their simulations, as they 
explicitly exclude stars from the ICL that got liberated from 
poorly resolved galaxies. Looking only at the three clusters 
they have simulated at high resolution, i.e. with a mass reso- 
lution comparable to our simulations, there is no significant 
discrepancy. For those three objects they find ICL fractions 
of 37%, 28%, and 41%, and no strong indication of a depen- 
dence on cluster mass. 



3.1.3 The radial distribution of BCG+ICL stars 

Fig. [7] compares the surface brightness profiles of the 
BCG+ICL component in thre e of our simu l ated groups to 
observational constraints from lZibetti et al.l l|2005h . For our 
simulated halos they were calculated for simulation outputs 
at redshift z = 0.21 of runs with and without A GN, scaled to 
the re ference redshift of 0.250 The profile from lZibetti et all 
(2005) is a mean surface brightness profile obtained by stack- 
ing SDSS cluster and group images, with a mean halo mass 
of M^ n « 7-8 x 10 13 M Q . Also shown is the profile they ob- 
tain without assuming a correction for incompletely masked 
satellite galaxies. In the figure, we compare with the two sim- 



7 The rp.25 - band of the photometric system defined 
IZibetti et al. I j2005l) is used for this comparison. 
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Figure 7. BCG+ICL surface brightness profiles of three groups 
simulated with and without AGN feedback are shown in the 
upper panel. For comparison observational constraints from 
Zibct ti et al ] l|2005l l obtained by stacking SDSS clusters and 
groups arc shown. Also shown is the mean surface brightness 
profile they find without correcting for incomplete masking of 
satellite galaxies (assuming a constant fraction of ~ 15% of un- 
masked galaxy light). For each simulated halo, the mass Mjqq 3 " is 
given in the figure's legend. The two groups shown i n red and blue 
have m asses very similar to the mean halo mass in lZibetti et al.l 
(2005). The lower pan e l show s the difference to the profile ob- 
tained in IZibetti et al.l l|2005l) . The vertical dotted line indicates 
the gravitational softening. 

to the mean mass of the observational sample. The mass of 
each simulated group is indicated in the figure's legend. Due 
to the removal of gas by AGN heating, the mass in the runs 
with AGN is typically slightly lower. In addition, we include 
results for a further, lower mass group. 

Without AGN, the surface brightness is typically too 
high in the central regions. In other words, the BCGs are 
overluminous in such simulations. When AGN are included, 
the central profile of one of the two groups, whose mass is 
close to the mean mass of the stacked clusters and groups, 
agrees very well with the observed central profile, indicating 
that a realistic BCG is formed in this simulation. The other 
massive group shown in the plot has a somewhat flatter 
central profile. However, given that there is object-to-object 
scatter, no perfect agreement for all halos close to the mean 
mass is expected. 

At radii larger than about 30 kpc, the surface brightness 
of the two more massive simulated groups exceeds the ob- 
served profile by roughly ~ 1 mag/arcsec 2 or somewhat less 
when AGN are included. This is consistent with our previ- 
ous findings that we have a very significant ICT component 
in our simulated halos, which as this comparison confirms 
appears to exceed the observed amount of ICL. As illus- 
trated further in the figure, we would need to compare to 
a simulated group of significantly lower mass to get a sur- 
face brightness profile that is close to the mean observed one 



within 200 kpc. However, in this case the simulated profile 
would then start to significantly underpredict the surface 
brightness at even larger radii. One should keep in mind 
though that for large radii this comparison may be quite 
sensitive to the accuracy of t he correction for in complete 
masking of satellite galaxies in IZibetti et ail (|2005T ). 

3.1.4 How do AGN affect the central cluster galaxies? 

Without a strong feedback mechanism, hydrodynamical sim- 
ulations of galaxy clusters suffer from excessive overcooling 
in cluster cores. Typically, the central galaxy is fed by a 
strong cooling flow and becomes too massive, too luminous 
and too blue compared to observations. One of the main 
aims of including AGN heating in cluster simulations was 
to offset cooling in cluster cores and produce more realis- 
tic central galaxies. Comparing the BCGs in runs with and 
without AGN feedback to observations should thus allow us 
to assess whether the employed AGN feedback model suc- 
cessfully regulates the central cooling flow. 

Fig. [5] shows the r-band luminosities of the BCGs in 
our simulations as a function of cluster mass. Results are 
shown for runs with and without AGN fee dback and are 
compa red to observational constraints from iPopesso et alj 
i|2007l ). Here Method 1 was used for finding the BCGs. The 
AGN heating successfully reduces the BCG luminosities for 
all halo masses and strongly improves the agreement with 
observations. For massive clusters the BCG luminosities ob- 
tained with AGN are still somewhat larger than observed. 
However, one should note that this depends somewhat on 
the details of the method used for making a distinction be- 
tween the BCG and ICL components in the simulations. 
When Method 3 or 4 is used, the resulting relation has much 
more scatter and the agreement with observations in runs 
with AGN is not as good as with Method 1 or 2. Note that 
the reduced BCG mass in our runs with AGN also reduces 
the gravitational lensing efficiency of our galaxy clusters (see 
iMead et alj|2010r ). 

3.2 When and where do the stellar populations of 
satellite galaxies, BCGs, and the ICL form? 

We now explore in what objects the stars have formed that 
are part of satellite galaxies, the BCG, and the ICL at red- 
shift z — 0. Using the IDs of the simulation particles, we 
can easily trace back every star particle in the simulation to 
the first snapshot after its formation. We then consider the 
SUBFIND (sub-)halo that the star particle is part of in this 
snapshot as its formation halo. Because we have saved sim- 
ulation snapshots very frequently, there is only a negligibly 
small number of particles that are not bound to any halo in 
the first snapshot after their creation. We can thus, in this 
way, assign formation halos to almost all star particles in 
the simulation. 

Fig. [9] shows the formation time of the star particles 
that are found in the BCG, the ICL, and the satellite mem- 
ber galaxies of a 10 14 /i _1 M© cluster at z = 0. Results are 
shown for simulations of this cluster with and without AGN 
feedback. Comparing the two runs clearly shows that AGN 
efficiently suppress star formation at low redshift. The ef- 
fect is particularly strong for the BCG where star formation 
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Figure 8. r-band BCG luminosity as a function of cluster mass 
Mjqq. Results are shown for runs with and without AGN feed- 
back. Method 1 (see Sect. 12, 21 was used for finding the BCGs. Ob- 
servational data_(meoUan BCG luminosities for six cluster mass 
bins) from [Popesso et al . (2007) are shown for comparison. The 
AGN feedback significantly reduces the BCG luminosities and 
improves agreement with observations. 



is almost completely shut off at redshifts z < 2. It is also 
interesting that in the run with AGN the BCG stars form 
on average earlier than the stars in the satellite galaxies, 
which nicely agrees with the observed cosmic "down-sizing" , 
where the mas s of galaxies hosti ng star formation decreases 
with time fe.g. lCowie et al 1996). At first, this observational 
finding is counterintuitive in a hierarchical structure forma- 
tion scenario . However, using sem i-analytic galaxy forma- 
tion models, |Pe Lucia et al.l 1 20061 ) have shown that AGN 
feedback allows reproducing this behaviour in a ACDM cos- 
mology where structure forms hierarchically. Our results in- 
dicate that the AGN feedback model in our cosmological 
hydrodynamical simulations operates in the same way. 

We also investigated in what objects stars form that at 
z = reside in a cluster's satellite galaxies, in the BCG, or 
the ICL. Fig. [10] shows the distribution of these stars accord- 
ing to the mass of the halo in which they have formed. Here, 
the halo mass is defined as the mass of the corresponding 
SUBFIND (sub-)halo. For the central galaxy of a halo this is 
roughly the same as A/200 i but for satellite galaxies it cor- 
responds to the smaller mass still bound within the gravi- 
tational tidal radius. The measured distributions are shown 
for a MfSo* = 10 14 /i _1 M Q cluster and for runs with and 
without AGN. 

Comparing the curves for the satellite galaxies, we can 
clearly see that AGN suppress star formation in halos larger 
than ~ 5 x 10 10 /i _1 Mq. Note that for smaller halos we do 
not expect any difference, as in our simulations black holes 
are only seeded in ha los that exceed this threshold mass (see 
iPuchwein et al.ll2008l ). This also means that in the runs with 
AGN the exact shape of the distribution at the low mass end 
will depend on the adopted threshold for the seeding. 

The stars that end up in the BCG at z — form in 
significantly more massive halos. This is not only due to 
star formation in the massive BCG itself, as for this clus- 
ter only ~ 10% and ~ 30% of the stars that reside in the 
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Figure 9. Distribution of the formation times of the stars in 
the satellite galaxies, the BCG, and the ICL of a 10 14 /i _1 Mq 
cluster. Results are shown for simulations with and without AGN 
feedback. Method 4 was used for making a distinction between 
BCG and ICL. The simulations were performed at zoom factor 3 
resolution (see Table [l}. 



BCG at z — have formed in the BCG's main progeni- 
tor in runs with and without AGN, respectively. Instead, 
also the stars that are acquired during mergers with the 
BCG tend to be formed in more massive halos. This is not 
unexpected, however, since the most massive galaxies are 
most likely to merge with the BCG due to their shorter 
dynamical friction timescale. We shall return to this point 
in Sect. 13.31 Without AGN, there is a significant fraction of 
BCG stars that form in very massive halos, i.e. in halos with 
masses larger than 10 13 /i _1 Mg. Indeed, many of these stars 
form in the cluster's main progenitor. On the other hand, 
in the run with AGN, star formation is almost completely 
shut off in the central galaxies of such massive halos, which 
again demonstrates that our model for AGN feedback can 
efficiently suppress strong cooling flows and excessive star 
formation in cluster cores. 

The curves for the ICL peak at roughly the same halo 
mass as for the satellite galaxies. However, the halo mass 
distribution is somewhat broader. There is also a second 
peak at very high halo mass, which is mostly due to stars 
formed in the cluster's main progenitor. We find a similar 
peak for most of our simulated clusters, with a position that 
depends on the cluster's mass. The occurrence of this second 
peak is rather surprising, especially since it is also found in 
runs with AGN, in which there is basically no star formation 
in the BCG once the cluster is as massive as required by the 
position of the peak. This means that these stars are not 
formed in the BCG but somewhere else in the cluster's main 
halo. 
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This is illustrated for the same cluster in more detail 
in Fig. 1111 which shows the distance of newly formed stars 
from the centres of their host halos. Stars that at z — 
are part of satellite galaxies or the BCG exclusively form 
in halo centres, roughly within 10/i _1 kpc. However, a sig- 
nificant fraction (~ 30%) of stars that are part of the ICL 
at z — form at much larger distances between 10/i -1 kpc 
and 1 /i -1 Mpc. We found a similar fraction also in our other 
clusters. 

We have investigated the nature of this 'intracluster star 
formation' in more detail, finding that it happens in small 
gas clouds that are distributed throughout the cluster, which 
can in fact be easily spotted in gas surface density maps. For 
most of these clouds, however, an associated halo can not be 
identified in maps of the dark matter particle distribution. 
We also checked whether these clouds form by thermal in- 
stability from the intracluster medium or whether they are 
stripped from infalling halos. We did this by tracing both the 
star-forming gas particles in these clouds, as well as all gas 
particles in the clusters main halo back to earlier times. This 
revealed that contrary to other gas particles, the vast ma- 
jority of star-forming gas particles fell into the cluster very 
recently. For example, less than 5% of the star-forming gas 
particles in these clouds in the progenitor of a 10 14 1i~ 1 Mq 
cluster at z = 1 are already part of this progenitor's main 
halo at z = 1.4, while this fraction is over 50% when look- 
ing at all main halo gas particles at the same cluster-centric 
radius. Also, at this earlier time most of the gas particles 
that later form intracluster stars can be unambiguously as- 
sociated with infalling halos. 

In other words, this means that most of these star- 
forming gas clouds consist of material that was stripped out 
of small infalling halos. Curiously, some of the gas clouds 
survived this stripping basically intact and continue to form 
stars, which then become part of the ICL. Note that there is 
observational evidence that some star formation does hap- 
pen in gas stripped from infalling galaxies (e.g. ISun et all 
Nevertheless, it is not entirely clear how realistic this 
mode of 'intracluster star formation' in our simulations is, 
and whether it is in part occurring due to numerical effects. 
For example, SPH si mulations are known to poorly resolve 
fluid instabilities fe.g. lAgertz et all 120071 1 that could act to 
disrupt such gas clouds after they are stripped from an in- 
falling halo. Also, thermal conduction in the hot intracluster 
medium is not included in our simulations and might play 
some role in reality, unless it is very efficiently suppressed 
by magnetic fields. 

If our simulations strongly overestimate the amount of 
such 'intracluster star formation', part of the discrepancy 
between our simulation results and th e amount of ICL in - 
ferred from observational studies like IZibetti et all (|2005l ) 
and lGonzalez et al.l (|2007l ) could be explained. In the most 
extreme case all star formation outside halo centres may be 
spurious, in which case we would expect about ~ 30% less 
ICL than reported above. On the other hand, it is also worth 
noting that if real clusters contain as many intracluster stars 
as predicted by our simulations, this would significantly re- 
lax the tension between constraints on the cosmic baryon 
fraction and baryo n fractions measured in galaxy clusters 
as, e.g., reported in lMcCarthv et al.l (|2007l ). Given the diffi- 
cult observational challenge to determine the amount of ICL 
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Figure 10. Distribution of the masses of the halos in which 
the stars in the satellite galaxies, the BCG, and the ICL of a 
10 14 ft. - 1 Mq cluster formed. Results are shown for simulations 
with and without AGN feedback. Method 4 was used for mak- 
ing a distinction between BCG and ICL. The simulations were 
performed at zoom factor 3 resolution (see Tabic [TJ. 



accurately, it certainly appears possible that the present ob- 
servations have still missed a significant amount of the ICL. 



3.3 In what objects do stars fall into the forming 
cluster? 

For each star particle that was not formed in the cluster's 
main progenitor we try to find the galaxy in which it fell 
into the forming cluster. We do this by looking for the most 
massive subhalo that a star particle was part of before be- 
coming part of the main cluster's FoF group. This allows a 
determination of the properties of these halos before tidal 
stripping in the cluster potential strongly affects them. 

Fig. [12] shows the distribution of the stars in satellite 
galaxy, BCG, and ICL components that were not formed in 
the cluster's main progenitor according to the time of their 
infall into the cluster. To reduce noise in the distributions we 
averaged them over four clusters, assigning the same weight 
to each of them. The simulations that were used for this 
figure included AGN feedback, but we note that simulations 
without AGN yield very similar results. The most obvious 
difference between the distributions is that almost all stars 
that end up in the BCG or ICL fall into the cluster before 
redshift z ~ 1, while the vast majority of stars that fall into 
the cluster later remain bound to their host galaxies. This 
suggests that galaxies falling into the cluster after redshift 
2 ~ 1 do usually not have enough time to sink towards the 
cluster centre by dynamical friction and are thus unlikely 
to merge with the BCG or get tidally disrupted close to 
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Figure 11. Distance of newly formed stars from the centres of 
the halos in which they were formed. Distributions are given for 
stars that at z = end up in satellite galaxies, the BCG, and 
the ICL of a 10 14 h~ 1 M Q cluster. Method 4 was used for making 
a distinction between BCG and ICL. Interestingly, a significant 
fraction of intraclustcr stars is not formed in halo centres. The 
increase of all distributions with distance for r values smaller 
than ~ 1 h~ 1 kpc is due to the increasing volume per radial bin. 



the centre. This further suggests that it takes a significant 
amount of time after the infall of a galaxy into a cluster 
until stars are efficiently stripped fro m it. Keeping t his in 
mind, it is not at all surprising that iMurante et alJ (|2007l ) 
find that most of the stripping happens at z < 1. 

We also checked the stellar masses of the galaxies in 
which the different components (satellite galaxy, BCG, and 
ICL stars) resided when they fell into the forming cluster. 
This is illustrated in Fig. 1131 The lower panel shows the dis- 
tribution of all infalling stars according to the stellar mass 
of the galaxy in which they fell into the cluster. The upper 
panel shows for each stellar mass bin of the infalling galaxy 
the fraction of stars that at z = reside in the cluster's satel- 
lite galaxies, the BCG, and the ICL. Again, the results were 
averaged over four clusters to reduce noise in the curves. 
For stellar masses below ~ 4 x 10 9 h~ 1 MQ, the fraction of 
stars that get stripped from an infalling galaxy and become 
intracluster stars strongly increases in our simulations. This 
is because these objects are poorly resolved, as also shown 
by our analysis of the convergence of the properties of the 
satellite galaxy population in Sect. 12.31 However, looking at 
the lower panel we see that only very few stars fall into the 
cluster in such small, poorly resolved galaxies. Their disrup- 
tion does therefore not significantly bias our predictions for 
the ICL. On the other hand, for more massive and well re- 
solved galaxies, we find that the fraction of stars that gets 
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Figure 12. Distribution (fraction per bin, each histogram is nor- 
malized to unity) of stars according to the time of their infall into 
the forming cluster. Results are shown for stars assigned to satel- 
lite galaxies, the BCG, and the ICL at redshift 2 = 0. Method 
4 was used for making a distinction between BCG and ICL. All 
curves were averaged over four clusters that were simulated with 
AGN. The same weight was used for each of them. 

liberated and joins the ICL increases with the galaxy's stel- 
lar mass. The most massive galaxies are also most likely to 
merge with the BCG, as can be seen by the strong increase 
of the BCG fraction. Accordingly, the fraction of stars that 
remains bound to an infalling galaxy decreases with stellar 
mass. This can be understood due to the shorter dynami- 
cal friction timescale for massive galaxies, while at the same 
time their tidal radius is larger. 

Furthermore, we investigated what fraction of intraclus- 
ter stars comes from completely dissolved galaxies. We did 
this by constructing a merger tree. For each SUBFIND (sub- 
)halo found in one snapshot we searched for the descendant 
(sub-)halo in the next snapshot that contains the largest 
number of the object's original star particles. In case this 
would yield the cluster's main halo as descendant, we also 
check which subhalo of the cluster contains the largest num- 
ber of the satellite galaxy's stars. If this subhalo still con- 
tains more than 10% of the original galaxy's star particles 
we consider it as the descendant, and not the cluster's main 
halo. This allows us to follow satellite galaxies even when 
they are strongly tidally stripped between two subsequent 
simulation snapshots. 

We find that the fraction of intracluster stars that comes 
from dissolved galaxies, i.e. galaxies that either merged with 
the BCG or were completely tidally disrupted, depends on 
halo mass. While we have performed this analysis only for 
four clusters, we find a systematic trend where this fraction 
drops from 75% for a 5 x 10 13 /i _1 M group to 22% for a 
3 x 10 14 Ii^Mq cluster. 

As discussed above, the mean fraction of stars that gets 
stripped from an infalling galaxy and becomes part of the 
ICL depends on the mass of the galaxy. On the other hand, 
this also means that the fraction of intracluster stars may 
be biased in simulations that do not reproduce the correct 
galaxy stellar mass function. Getting the latter right in cos- 
mological hydrodynamical simulations is, however, a long- 
standing problem which has not yet been solved satisfac- 
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Figure 13. Distribution (fraction per bin) of stars according to 
the stellar mass of the galaxy in which they fell into the forming 
cluster (lower panel). The upper panel shows for each stellar mass 
bin of the infalling galaxies the fraction of stars that end up in 
satellite galaxies, the BCG, and the ICL at redshift z = 0. Method 
4 was used for making a distinction between BCG and ICL. All 
curves were averaged over four clusters that were simulated with 
AGN. The same weight was used for each of them. 



torily. Looking at the upper panel of Fig. 1131 we see that 
we might overestimate the ICL fraction if we significantly 
overpredict the number of very low mass or very massive 
galaxies. The former does not seem to be a problem as in- 
dicated by the distribution in the figure's lower panel. The 
abundance of massive galaxies, however, is typically overpre- 
dicted in cosmological hyd rodynamical simulations (see e.g. 
lOppenheimer et al.l (|2009l )). This surely is the case in our 
simulations without AGN, but we can not be fully sure yet 
whether this is completely resolved when the AGN feedback 
model is included. Assessing this reliably requires a high res- 
olution simulation of a large cosmological box, which has not 
yet been done due to the large computational cost. 

Nevertheless, from exploring the cluster galaxy mass 
functions of our simulated clusters, we can be confident 
that the problem of overpredicting massive galaxies is rather 
strongly alleviated in runs with AGN, as also suggested by 
the much lower BCG luminosities shown in Fig. [8] However, 
the ICL fractions in our simulations with AGN feedback are 
even slightly higher than in the runs without it. This sug- 
gests that overpredicting massive galaxies is not the main 
reason why our ICL fractions are relatively high compared 
to the values typically inferred from observations. This also 
highlights that it is really far from obvious how a large ICL 
fraction can be avoided in the simulations of the formation 
of clusters in the ACDM cosmology. 



3.4 The fate of infalling galaxies 

Fig.[l4]illustrates the infall of satellite galaxies and the origin 
of intracluster stars for a 10 14 h~ x MQ cluster. The upper left 
panel shows all star particles residing in the cluster's satellite 
galaxies, and the BCG and ICL components at z — 0. Here, 
Method 4 was used for making a distinction between BCG 
and ICL, and the figure illustrates the assignment of the 
star particles to the individual components. The upper right 
panel of Fig. [14] shows the same cluster at z = 1. All star 
particle are still coloured according to the component which 
they belong to at z = 0, allowing some insights into the 
assembly history of the cluster galaxy population and the 
origin of intracluster stars. 

Looking at some of the infalling galaxies, it can be eas- 
ily seen that already at z — 1 stars that later become part of 
the ICL are preferentially found in the outskirts of infalling 
galaxies, and are therefore comparatively weakly bound to 
them. Furthermore, the figure shows that only the most mas- 
sive infalling object (marked as A) contributes a significant 
number of stars to the BCG, as expected from our analy- 
sis in the previous section. Following the evolution of clus- 
ters in this way also nicely shows that galaxy mergers often 
create a loosely bound component of stars which is sub- 
sequently stripped when the merger remnant falls into the 
cluster. Unfortunately, it is a bit hard to il lustrate this in 
a still image, but it confirms the finding of iMurante et all 
(2007) that mergers in the assembly history of the BCG and 
of other massive cluster galaxies are of critical importance 
for the formation of the ICL. 

We now consider the fate of some infalling galaxies in 
more detail. For this purpose we have selected two infalling 
objects marked as A and B in the upper right panel of Fig. 1141 
and followed them to redshift z = 0. The distribution of 
their stars at z — is shown in the lower panels. The core 
of group-sized object A has basically merged with the BCG 
at this point. However, even after several passages of the 
cluster core one can still see significant tidal features both 
in the distribution of stars that at z — are assigned to 
the BCG and the ICL. The infalling galaxy B, on the other 
hand, remains largely intact and becomes one of the cluster's 
satellite galaxies at z = 0. Those stars that were stripped 
from it form a large tidal stream extending almost over the 
whole cluster. These two examples illustrate that the ICL 
really consists of many individual tidal streams and shell- 
like features. Only their superposition looks like a smooth 
distribution of intracluster stars. 

As discussed above, stars that later become part of 
the ICL are preferentially found in the outskirts of infalling 
galaxies (see also the upper right panel of Fig. I14[l . It also 
seems plausible that the most weakly bound stars are the 
ones most likely to be stripped from an infalling galaxy and 
to become intracluster stars. We have checked this explicitly 
by calculating the binding energy of stars in such infalling 
galaxies, and Fig. [TS] illustrates the results of our analysis. 
We show the fraction of stars that at z = reside in the 
cluster's satellite galaxies, the BCG, and the ICL as a func- 
tion of their binding energy in objects A and B at z = 1. 
As expected, the fraction of stars that become part of the 
ICL strongly decreases with increasing binding energy. For 
object A, the most bound stars end up in the cluster's cen- 
tral galaxy, while for object B, the most bound stars remain 
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Figure 14. Distribution of stars in a 10 14 h~ 1 Mq cluster simulated with AGN feedback. The upper left panel shows the stars residing 
in the cluster's satellite galaxies (blue), BCG (green), and ICL (red) at redshift z = 0. Grey dots indicate stars that are not part of 
the clusters FoF group at z = 0. Method 4 was used for making a distinction between BCG and ICL. The upper right panel shows the 
cluster's progenitors at z = 1. All star particles are coloured as before, i.e. according to the component in which they end up at z = 0. 
The circles mark two infalling objects A and B. In the lower panels the positions of the stars that are part of A (lower left) and B (lower 
right) at z = 1 are shown at z = 0. Again they are coloured according to the component to which they are assigned at z = 0. All stars 
that are not part of A or B are show in grey in the lower panels. 



bound together in a satellite galaxy. Overall, about half of 
the stars in A end up in the BCG, while the other half be- 
comes part of the ICL. For object B, only about ~ 10% of 
the stars are stripped and become intracluster stars. 

We have shown in Fig. [13] that the fraction of stars 
that gets stripped from an infalling galaxy depends on the 
galaxy's mass. However, it also depends on the orbit of the 



galaxy. This is illustrated in Fig. 1161 which shows the frac- 
tion of stars that at z = are still bound to a satellite galaxy 
as a function of the minimum distance of the galaxy's or- 
bit from the cluster centre. The fractions are shown for all 
galaxies that fell into the progenitor of a 3 x 10 14 h~ 1 M(? ) 
cluster between redshifts z — 1.5 and z = 1. In order to de- 
termine the minimum distance, the orbit of the galaxy and 



© 0000 RAS, MNRAS 000, 000-000 



Intracluster stars in simulations with AGN feedback 15 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

binding energy in units of the most bound particle 

Figure 15. Distribution of binding energies in their host halos 
at z = 1 of the stars that end up the cluster's satellite galaxies, 
BCG, and ICL at z = 0. Results are shown for the stars that are 
part of the objects A (upper panel) and B (lower panel) indicated 
in Fig. 1141 Method 4 was used for making a distinction between 
BCG and ICL. 



of the cluster centre were interpolated using all simulation 
snapshots between the time of the galaxy's infall and z — 0, 
using a cubic interpolation between each pair of successive 
snapshots based on the halo positions and velocities deter- 
mined by SUBFIND. The results are shown for three different 
stellar mass ranges of the infalling galaxy. We see that only 
few stars are stripped from galaxies that never come close to 
the cluster centre. On the other hand, galaxies that closely 
approach the cluster centre are often substantially stripped, 
completely tidally disrupted, or merge with the BCG. This 
confirms the expectation that most of the tidal interactions 
that liberate intracluster stars happen close to the cluster 
centre and the BCG. 
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Figure 16. Fraction of stars that remain bound to a satellite 
galaxy as a function of the minimum distance of its orbit from 
the cluster centre. All galaxies that fell into the main progenitor 
of a 3 X 10 14 k _1 Mg cluster between redshift z = 1.5 and 2 = 1 
are plotted. The minimum distance was found by interpolating 
the galaxy's orbit based on all simulation snapshots between the 
time of the infall and z = 0. Different symbols are used for low 
mass galaxies (stellar mass before infall into the cluster M, in the 
range 5 X 10 9 h _1 M < M, < 5 X 10 10 /i _1 M Q ), intermediate 
mass galaxies (5 X 10 10 h^ 1 M Q ^ M* < 5 X 10 11 h~ 1 M e ) and 
high mass galaxies (M„ ^ 5 X 10 11 /i -1 Mq). Galaxies with a 
fraction of have either merged with the BCG or have been 
completely tidally disrupted. 



4 SUMMARY AND CONCLUSIONS 

The intracluster light in clusters of galaxies represents an 
interesting and significant component of their total stellar 
mass. While the observational constraints on the ICL are 
still uncertain, some studies have reported intracluster star 
fractions of up to ~ 50%. Even if the true fraction is sig- 
nificantly lower, it is therefore clear that the ICL can not 
be neglected in the baryon and stellar mass budgets of clus- 
ters. Furthermore, the radial profile of the ICL and the mass 
contained in it may pose interesting constraints on galaxy 
cluster formation models. Yet, comparatively little theoret- 
ical work has been carried out thus far on the formation of 
the ICL. In fact, most semi-analytic models of galaxy for- 
mation have hitherto ignored this component entirely. 

In this work we have therefore studied the origin of in- 
tracluster stars in a set of high-resolution hydrodynamical 
simulations of the formation of clusters of galaxies, embed- 
ded in their appropriate cosmological setting. Our sample 
of resimulations has been randomly drawn from the Mil- 
lennium simulation, with the only selection criterion being 
to provide a wide coverage of group and cluster masses, of 
roughly two orders of magnitude in halo mass. Thanks to the 
very high mass and force resolution of our resimulations, our 
simulations provide a powerful and representative sample of 



the whole cluster population. Another very timely aspect 
of our simulations is that we not only account for hydro- 
dynamics, radiative cooling, heating by a UV background, 
star formation and supernovae feedback, but also incorpo- 
rate a state-of-the-art model for the growth of supermassive 
black holes and for feedback from AGN. Because we simu- 
lated each cluster both with and without AGN physics, this 
allows us to pinpoint the impact of AGN heating on the 
cluster galaxy populations, and in particular on the ICL. 

Our results clearly confirm the importance of AGN 
feedback in galaxy clusters. AGN lead to a reduction of 
the amount of stars in our clusters and groups by about 
one third, roughly independent of cluster mass. Especially 
the stellar masses and luminosities of BCGs are greatly re- 
duced by AGN feedback, and their stellar populations be- 
come much older. As a consequence, the BCGs are in much 
better agreement with observational constraints. Further- 
more, in poor clusters and groups, the total baryon fractions 
within r§oo become significantly lower when AGN heating is 
included. 

The primary focus of our analysis has been on the 
amount and the origin of the ICL component in our sim- 
ulated groups and clusters. In order to allow a meaningful 
comparison of our simulation results with observations of 
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clusters and groups, it was necessary to find a robust way to 
assign star particles in the simulations either to one of the 
cluster's satellite galaxies, to the BCG, or to the ICL. Espe- 
cially the distinction between the latter two components is 
not without ambiguities. We therefore developed and tested 
several different methods for distinguishing between the dif- 
ferent components. This helps to assess the robustness of 
the measurements and allows the validation of measurement 
techniques that are close to the methods usually employed 
in observational studies. 

Our investigation of the properties of the satellite 
galaxy populations, the BCGs and the ICL components in 
our simulations yielded the following main results: 

• We find a very significant fraction of ~ 40% (without 
AGN) or ~ 45% (with AGN) of intracluster stars within 
rloo in our simulated clusters and groups. These values are 
robust with respect to numerical resolution and integration 
accuracy, and are almost independent of halo mass. They 
are, however, larger than those typically inferred from ob- 
servations. 

• For all the different methods we tested for making a 
distinction between BCG and ICL, we find that the vast 
majority of stars in the main halo of our simulated clusters 
are part of the ICL, rather than the BCG. As a result of this 
dominance, our ICL fractions also do not strongly depend 
on the exact method used for making the distinction. 

• The intracluster stars form on average earlier than the 
stars residing in cluster satellite galaxies. Stars ending up in 
the BCG at z = 0, however, typically form at even earlier 
times, especially when AGN are included in the simulations, 
which prevent excessive star formation in BCGs at low red- 
shift. 

• The fraction of stars that is stripped from infalling 
galaxies increases with galaxy stellar mass, which can be 
understood as a consequence of their larger ratio of size to 
tidal radius. The orbit of an infalling galaxy also plays an 
important role. Galaxies that closely approach the cluster 
centre typically loose the largest fraction of stars. 

• We find that most intracluster stars are stripped from 
galaxies that fell into the cluster early, roughly before z ~ 1. 
Most likely, galaxies falling in later just do not have enough 
time to sink towards the cluster centre by dynamical friction 
and to get tidally disrupted there. 

• The stars stripped from individual galaxies are not 
smoothly distributed within the cluster even after several 
orbits of their former host galaxy, instead they form streams 
and other tidal features. The ICL is therefore not smooth 
but a superposition of many such tidal features. 

• In our simulations, a significant fraction of intracluster 
stars (up to ~ 30%) do not form in galaxies but in cold gas 
clouds stripped from infalling substructures. These clouds 
remain intact after being stripped and give rise to 'intra- 
cluster star formation'. 

• When including AGN, the total stellar masses in our 
simulated galaxy groups are in good agreement with obser- 
vations. However, the simulations do not reproduce the steep 
observed decline of the stellar mass fraction with halo mass, 
so that our massive simulated clusters contain more stars 
than suggested by observations. 

It will be very interesting to see whether improved ob- 
servational determinations of the ICL fraction in clusters 



confirm a lower value than found in our simulations. This 
would represent a non-trivial challenge for future simulation 
models, as our work has shown that the simulated ICL frac- 
tion is remarkably robust, not only with respect to integra- 
tion parameters and numerical resolution, but also to quite 
drastic changes in the modelling of star formation. This also 
means that it is not obvious at all how a cluster simulation 
could be obtained that reproduces the galaxy population as 
well as in our best runs but at the same time yields a much 
lower ICL component. 

One interesting effect we found in our simulations is 
that a sizable fraction of our intracluster stars actually forms 
'in situ', in gas clouds that were stripped out of the dark 
matter halos of infalling galaxies but that are not bound 
to a dark matter satellite any more. It is not entirely clear 
whether this mode of 'intracluster star formation' is signifi- 
cantly affected by numerical inaccuracies in the simulations. 
For example, fluid instabilities might be able to disrupt these 
clouds and suppress the star formation in them if they are 
better resolved than possible with SPH. If this were indeed 
the case, this could explain some part of the discrepancy. 
On the other hand, a high ICL fraction would actually help 
to relax the current tensions between constraints on the cos- 
mic baryon fraction and the baryon fractions measured in 
galaxy clusters. 
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